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Experiments with two-dimensional (2D) plasma crystals are usually carried out in rf plasma sheaths, where the 
interparticle interactions are modified due to the presence of plasma wakes. The wake- mediated interactions 
result in the coupling between wave modes in 2D crystals, which can trigger the mode-coupling instability 
and cause melting. The theory predicts a number of distinct fingerprints to be observed upon the instability 
onset, such as the emergence of a new hybrid mode, a critical angular dependence, a mixed polarization, 
and distinct thresholds. In this paper we summarize these key features and provide their detailed discussion, 
analyze the critical dependence on experimental parameters, and highlight the outstanding issues. 



PACS numbers: 52.27.Lw, 52.27.Gr 

I. INTRODUCTION 

Stability of two-dimensional (2D) plasma crystals is a 
fundamental problem of complex (dusty) plasmas. Such 
crystals - which are monolayers of hexagonally ordered 
monodisperse microparticles - can be (routinely) created 
in a rf plasmai*^: Particles get negatively charged in a 
plasma and therefore the electric force exerted on them in 
the sheath above a horizontal electrode can compensate 
for gravity, thus providing a stable levitation. There are 
several mechanisms operating in complex plasmas that 
can result in the melting of 2D crystals. These mech- 
anisms can generally be divided into two categories - 
generic and plasma- specific. 

Generic mechanisms of melting are those operating in 
any (classical) system with a given pair interaction be- 
tween particles, provided the interaction can be described 
by a Hamiltonian. The melting in 2D systems can ei- 
ther be a two-step process (which involves consecutive 
unbinding of weakly interacting dislocation and discli- 
nation pairs, respectively, with the intermediate hexatic 
phase)i2r£ or a one-step process (where the hexatic phase 
is preempted by the formation of dislocation chains)^. 
These generic melting mechanisms can operate in very 
different 2D systems including complex plasmas^—. 

Plasma-specific melting mechanisms can only operate 
in complex plasmas. Such mechanisms are associated 
with the energy exchange between microparticles and 
ambient plasma and can be considered as a result of the 
system openness. For instance, 2D plasma crystals can 
be strongly perturbed by single particles moving above 
or below the monolaye r or they can melt due to fiuc- 
tuations of particle charge^i^— . 

The most universal among the plasma-specific mech- 
anisms is that associated with the wake-mediated inter- 
action between microparticles: In the presence of strong 
plasma flow the screening cloud around each charged par- 
ticle becomes highly asymmetricJ^.^^. These clouds are 
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usually referred to as "plasma wakes"— — and play the 
role of a "third body" in the interparticle interaction, 
making it nonreciprocal^^. Under certain conditions, this 
makes the system non-Hamiltonian and provides effective 
conversion of the energy of flowing ions into the kinetic 
energy of microparticle^^ i^^i^^ . 

The wake-induced mechanism of melting of crystalline 
monolayers was discovered theoretically a decade ago by 
Ivlev and Morflll^^. Based on a simple model of a par- 
ticle chain, it was shown that the longitudinal in-plane 
and transverse out-of-plane dust-lattice (DL) wave modes 
are no longer independent - they are coupled due to the 
wake- mediated interactions. When the modes intersect 
they become modifled by the coupling and form a hybrid 
mode in a narrow vicinity of the crossing. This can trig- 
ger the mode- coupling instability which causes the melt- 
ing. This melting mechanism had received strong conflr- 
mation later on, when the instability threshold predicted 
by the theory was compared^^ with the experimental ob- 
servations by U. Konopka and numerical simulations by 
G. Joyce. However, the direct comparison of theory and 
experiment became possible only recently, after an ex- 
perimental method of measuring the out-of-plane mode 
was developed^^ and therefore the essential flngerprint of 
the mode-coupling instability - the hybrid mode - be- 
came observable. In recent experiment by Couedel et 
al^ an implementation of this method unambiguously 
demonstrated that the melting indeed occurs due to the 
resonance coupling between the longitudinal in-plane and 
transverse out-of-plane modes. The variation of the wave 
modes with the experimental conditions, including the 
emergence of the hybrid branch, revealed exceptionally 
good agreement with the theory of mode-coupling insta- 
bility generalized for 2D case by Zhdanov et alS^. 

The theory of mode-coupling instability^iiHi^ predicts 
a number of distinct flngerprints to be observed upon the 
instability onset: Along with the emergence of hybrid 
mode mentioned above, these are a critical angular de- 
pendence - the hybrid mode flrst appears only for wave 
vectors oriented along one of the principal lattice axes, 
a mixed polarization - the two DL modes that form the 
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hybrid mode are no longer purely longitudinal and trans- 
verse close to the merging point, and distinct thresholds 
- the instability sets in only when (i) the particle number 
density in the monolayer is high enough or /and vertical 
confinement eigenfrequency is low enough (so that the 
two DL modes can cross and form the hybrid mode) and 
(ii) the gas pressure is low enough (so that the growth 
rate of the hybrid mode exceeds the damping rate). 

All these fingerprints have been mentioned in our pre- 
vious theoretical publication^i^i^, some of them were 
also illustrated in the followup experimental paper—. 
However, their discussion was apparently too concise to 
address all these important properties in necessary detail. 
The need for in-depth discussion of the wake-mediated 
mode coupling became evident now, when new publica- 
tions have appeared where some essential properties of 
the instability were misinterpreted (see, e.g., recent ex- 
perimental paper^^ and the subsequent Erratum^^ by Liu 
et al.). Therefore, in this paper we summarize the key 
features of the mode couphng and provide a detailed dis- 
cussion, analyze the critical dependence on experimental 
parameters, and highlight the outstanding issues. 

II. SUMMARY OF MODE-COUPLING THEORY 

In 2D plasma crystals, two in-plane wave modes can 
be sustained (here we naturally leave aside polydisperse 
mixtures). Both modes have an acoustic dispersion, one 
of them is longitudinal and the other is transverse. Since 
the strength of the vertical confinement in such systems 
is finite, there is a third fundamental wave mode which 
has an optical dispersion and is associated with the out- 
of-plane oscillation a^-'^i^^i^^ — . Theory predicts that all 
three modes depend critically on the parameters of the 
plasma wakes. Yet the role of the wakes is not only the 
modification of the modes themselves. What is far more 
important is that the modes are coupled to each other 
due to the wake- mediated interactions between particles. 
It is worth noting that this coupling is linear and there- 
fore does not depend on the wave amplitude. Below we 
identify the conditions when the wake-mediated coupling 
becomes crucial for waves in 2D crystals. 

A. DL modes for nonreciprocal interactions 

In general, linear dispersion relations c<j(k) are deter- 
mined by eigenvalues (eigenmodes) of a dynamical matrix 
D. The latter is derived by considering small perturba- 
tions of individual particles (with respect to a stable con- 
figuration, e.g., a hexagonal monolayer for 2D crystals) 
of the form oc exp(— zcjt + ik • r). Elements of D are de- 
termined by properties of interparticle interactions and 
are functions of k. If the interactions are reciprocal then 
D is Hermit ian, and therefore the eigenvalues are always 
real (i.e., the modes are stable) and the eigenvectors are 
orthogonal. The situation changes if the interactions are 




FIG. 1. (a) Elementary hexagonal lattice cell with the frame 
of reference, the lattice constant is A. (b) The reciprocal 
lattice in k-space, the basis vectors of the lattice are bi,2 = 
27rA~"'^(l, i-^)- Due to the lattice symmetry, it is sufficient 
to consider the wave vectors kat0°<^<30° and from 
within the first Brillouin zone (gray region enclosed by dashed 
lines), so that |k|A < fyr for 6' = 0° and |k|A < ^tt for 
= 30°. 

nonreciprocal - the eigenvalues in this case can become 
complex. 

Before specifying a particular functional form of in- 
terparticle interactions in 2D crystals, let us generally 
assume the non-reciprocity and discuss common proper- 
ties of the DL waves in this case. The dynamical matrix 
has the following form: 

/ ah - /3 27 iay \ 

D = 27 ah + /3 id, . (1) 

\ iay ia^ nl^^^ - 2av / 

The elements ah(k), /3(k), and 7(k) determine the dis- 
persion of two in-plane (horizontal) modes and av(k) 
characterizes the out-of-plane (vertical) mode. The el- 
ements o-x,y(k) emerge due to non-reciprocity and make 
D non-Hermitian. The matrix is calculated for the ref- 
erence frame shown in Fig. [T]^a) assuming that vertically 
the particles are confined in a (parabolic) potential well 
characterized by the eigenfrequency l^conf = 27r/conf- To 
derive the dispersion relations a;(k), we add a friction 
force miyr in equations of motion (with u being the damp- 
ing rate) which yields: 

det[D - uj{uj ^ iu)\] = 0, (2) 

where I is the unit matrix. Thus, uj{uj -\- iv) = are the 
eigenvalues of D, i.e., the DL wave modes. 

For reciprocal interactions cr^^y =0 and then Eqs ([T]) 
and (|2|) yield 

{n^-^l^^){n^-^l^){^''-^l) = Q. (3) 

There are three independent and real DL modes: Two 
acoustic in-plane modes Q?{k) = ah ± + 47^ = 
^h|| (where || and ± indicate the longitudinal and 

transverse polarizations and correspond to the plus and 
minus signs, respectively) and an optical out-of-plane 



mode Q'^{k) = Q^^^^ - 2a^ = nl{k). When solved for 
cc;(k), each mode yields a couple of conjugate branches. 

When interactions are nonreciprocal the modes are 
modified - they are described by Eq. ([3]) with nonzero 
right-hand side (r.h.s.) which is proportional to + 
The DL modes become coupled with each other and 
comple:^. 

Now let us specify interparticle interactions. For the 
direct interaction between charged grains we naturally 
adopt the Yukawa potential which is characterized by 
(negative) particle charge Q and the screening length 
A. For the wake we employ the commonly used simple 
model^ of a point-like (positive) charge q located at the 
distance 5 downstream each particle and also interacting 
with the neighboring particles via the Yukawa potential. 
The elements of the dynamical matrix ([1]) for such inter- 
actions are calculated in Appendix [Aj Eqs ()Aip and ()A2p . 
For brevity, theoretical results in this paper are written 
in the dimensionless form - the frequency is normalized 
by the DL frequency scale, 




and the wave vector is normalized by the lattice constant 
A: ft/n^L ^ and kA k. 

B. Types of mode coupling 

After the normalization, the DL modes depend on 
three dimensionless parameters (along with the normal- 
ized wave vector): The screening parameter = A/A 
as well as the relative wake charge q = \q/Q\ and dis- 
tance S = 6/ A. In experiments, is typically about 
unity whereas both wake parameters are usually small. 
Elements (Jx^y [Eq. ()A2p ] are proportional to the prod- 
uct qS and hence the coupling between the DL modes is 
very small [ex (qS)'^]. It becomes important only if dif- 
ferent DL modes intersect (or they are very close to each 
other); otherwise the coupling can be neglected and then 
the modes can be treated independently, as described by 
Eq. (|3]). Let us elaborate on this point. 

Behavior of two coupled modes, l]i^2(k), in a close 
proximity of their intersection can generally be described 
by the following equation: 

[1] - 1^0 - Ui . (k - ko)] [n-no-U2-{k- ko)] = e. 

(4) 

Here Ui,2 = 91^i,2/Sk are the group velocities of the 
two modes calculated at a point ko of the intersection 
line 1^0 (k) (the latter is determined by the shape of the 
crossing modes) and e characterizes the coupling. From 
Eq. (|4j) one readily concludes that the effect of the cou- 
pling critically depends on the sign of e. This is illus- 
trated in Fig. [2) When e < the crossing modes merge 
and form a complex hybrid mode in a narrow range 
I k - ko I (X v/e, with Re l^hyb r^o + | (Ui + U2 ) • (k - ko ) 
and Im l^hyb ^ ^^/^' Thus, one of the branches of the 
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FIG. 2. Modification of coupled modes [Eq. Q] in the vicin- 
ity of their crossing. The sketch depicts the modes in the 
(Q, |k|)-plane (a) in the absence of coupling, (b) with nega- 
tive coupling (e < 0), and (c) with positive coupling (e > 0). 
Negative coupling results in the formation of a hybrid mode 
which can be unstable (see text); positive coupling causes the 
mode reconnect ion, while the modes remain stable and form 
a "forbidden band". 



hybrid mode has a positive imaginary part and hence can 
be unstable. When e > the coupling is accompanied by 
reconnection (formation of a "forbidden band") and the 
modes remain stable (real). Furthermore, the sign of the 
product Ui • U2 in the unstable case determines the type 
of the instability, which is absolute for Ui • U2 < (the 
case illustrated in Fig. [2]) and convective otherwise. 

Now we can consider the DL modes described by Eqs 
([1]) and ([2|), with the dispersion and coupling elements 
from Eqs (|Aip and (|A2p . The detailed analysis shows 
that the intersection of the transverse in-plane mode 
^h±(k) with the longitudinal in-plane, l]h||(k), as well 
as with the out-of-plane mode l^v(k) is characterized by 
positive coupling and therefore is accompanied by recon- 
nection illustrated in Fig. [Sh, so that the modes remain 
stable. Moreover, the coupling for these two pairs of 
modes exactly disappears (Fig. [2^) if k is parallel to one 
of the principal lattice axes (i.e., when ^ = 0° or 30° in 
Fig.H^). 

In contrast, for the remaining pair of l^v(k) and l^hy (k) 
the coupling is always negative. As illustrated in Fig. [2Jd, 
this results in the formation of the hybrid mode with the 
imaginary part Im l^hyb ^ (by the order of magni- 
tude). If damping is low enough, this triggers the mode- 
coupling instability. 

Thus, only the intersection between the out-of-plane 
and longitudinal in-plane modes is critical for stability of 
2D plasma crystals^. 



C. Mode-coupling theory for ID model 

It is instructive to present a brief summary of the 
mode-coupling theory which was originally developed 
for a ID string assuming the nearest-neighbor (NN) 
interactional^. As it follows from the subsequent 
analysis (see Sec. IIIip . this approach yields very sim- 
ple expressions which can be conveniently implemented 
for practical purposes (e.g., to calculate the instability 
thresholds, see Sec. IIIID"]) . 

In the framework of ID NN modeP^, the coupling be- 
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FIG. 3. Fluctuation spectra for "shallow" mode intersection (molecular dynamics simulations). Shown are the DL modes 
(positive branches) for k at (a) ^ = 0° and (b) — 30°, and (c) spectrum in the k-plane integrated over frequency (in the 
range between 10 Hz < / < 14 Hz). Simulations are for N — 16384 particles of a mass m = 6.1 x 10"^"^ kg and charge 
Q — —18500 e, the screening length is A = 600 /xm, the damping rate is v — 0.87 s~^. The horizontal and vertical confinements 
have eigenfrequencies f^\^ — 0.11 Hz and /^o^f — 14.03 Hz, respectively. The lattice constant in the center is A ~ 612 /xm 
{hi — A/ A :^ 1.02). The wake parameters are q — —0.2(3 and 5 — 0.3 A. The dashed lines show the border of the first Brillouin 
zones, solid lines in (a) and (b) are theoretical curves. In (c) (as well as in other figures depicting spectra in the k plane) the 
calculations for the upper right quadrant are mirrored in other three (shown in lighter colors). 



tween horizontal ("longitudinal in-plane") and vertical 
( "out-of-plane" ) DL modes, 

^lik) = 4(1 - q) (/^-i + 2/^-2 + 2/^-3) e-'" sin^ \k, 
^Uk) = ^c'onf - 4(1 - q) (^-' + ^-') sin^ ^k, 

(5) 

is described by the following (exact) equation: 

{n^ -nl){n^ -nl)^a^ = 0. (6) 

The range of k is limited by the first Brillouin zone \k\ < 
7T (we use the same normalization as above), and the ID 
coupling coefficient is 

a{k) = 2qS {k'^ + 3/^"^ + 3/^"^) e"'" sin k. (7) 

All characteristics of the hybrid mode (e.g., width, 
growth rate, etc.) can be easily derived from Eqs (j5j)- 
([7]), see Appendix [Bl 

These results are obtained in the limit of "small" S 
(for general case, see Ref^^). In Sec. IIII Dl thev will be 
compared with the exact 2D results derived in this limit. 



D. Additional remarks 

To conclude this Section, we should make two spe- 
cific remarks which might be useful when considering DL 
waves and their coupling: 

(i) The dynamical matrix D is solely characterized by 
interparticle interactions, i.e., does not depend on the 
damping rate v. Thus, the wave modes ^7^(k) - which 
are the eigenvalues of D - do not depend on v either. As 
regards the dispersion relations ct;(k), their dependence 
on u is determined by the solutions of the quadratic equa- 
tion uj{uj -\- iu) — = 0. Therefore the effect of friction 
on waves is straightforward: For weakly damped waves 
(when \Q\ ^ z^, which is typical for experiments), one 
readily obtains co'(k) c:^ l](k) — |iz/. From the practical 
point of view this implies that one can analyze undamped 
dispersion relations (i.e., put uj = Q) and afterwards sim- 
ply add —^v to the resulting imaginary part^, which si- 
multaneously identifies the instability threshold for the 
hybrid mode"^^ (see Sec. IIII D 2|) . 

(ii) For the analysis of wave modes in a crystal it is 
sufficient to consider the wave vectors from within the 
first Brillouin zone&^ which is shown in Fig. [Ud. This 
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zone is nothing but the Wigner-Seitz ceh of the reciprocal 
lattice formed by the basis vectors bi and b2. Hence, the 
wave vectors k and = k + G which are different by a 
linear combination of the basis vectors (G = mbi -\-nh2) 
are equivalent for wave modes, i.e., l](k + G) = l](k). 

One has to remember this fact when analyzing dis- 
persion relations at large |k|. For instance, recently Liu 
et alS^ identified one of the "hot spots" seen in their 
fluctuation spectra (in their Fig. 4c) as a new hybrid 
mode. This spot is centered at |k|A 6.3(::^ 27r), 6> = 0° 
and therefore is located outside the first Brillouin zone; 
after the mapping with, say, G = b2, its position is 
|k|A ^TT, e = 90° (equivalent to 30°, see Fig. [T]d; 
note also the change in the polarization, from the trans- 
verse to longitudinal). This naturally coincides with the 
position of the "regular" hybrid mode at the border of 
the first zone, whereas the "original" spot is merely its 
image outside it (as illustrated in FigsjSh andHJ)). 

III. FINGERPRINTS OF MODE COUPLING 

In this Section we present detailed discussion and 
analysis of distinct fingerprints characterizing the wake- 
mediated coupling of DL modes in 2D crystals. Based on 
the results of Sec. lIIBl we naturally focus only on the cou- 
pling between the out-of-plane and longitudinal in-plane 
modes, which causes the formation of the hybrid mode 
and can trigger the mode-coupling instability. 

To illustrate key features of the mode coupling, the 
theoretical calculations were combined with the results 
of molecular dynamics (MD) simulations (see descrip- 
tion in Appendix [C]) . The DL wave modes were obtained 
from the simulations by plotting the so-called "fluctua- 
tion spectra" - intensity distributions of thermal velocity 
fluctuations in the (a;,k) space. 

A. Hybrid mode 

Figure [3] represents a characteristic example of fluctua- 
tion spectra measured at the onset of the mode-coupling 
instability. It was obtained from MD simulations per- 
formed for the conditions close to those of experiment by 
Couedel et al.^^ (experiment II in Table H]). In this case 
the out-of-plane and longitudinal in-plane modes are just 
barely touching. Such "shallow" intersection is never- 
theless sufficient to form a distinct hybrid mode and to 
trigger the instability. 

Figures [3^ and b depict the DL modes in the (cj, |k|) 
plane (positive branches). They are obtained for two di- 
rections of k at ^ = 0° and = 30°, respectively (corre- 
sponding to the principal lattice axes, see Fig. [TJi). Fig- 
ure [3t shows the intensity distribution in the k-plane (in- 
tegrated over frequency) , which provides a "top view" on 
the fluctuation spectra (dotted line shows the direction of 
= 30°). The hybrid mode can be identified in Figs[3]3 
and c as "hot spots" located near the border of the first 
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FIG. 4. (a) Contour plot of the growth rate Im C(;hyb(k) calcu- 
lated from theory and (b) fluctuation spectra in the k-plane 
obtained from simulations (integrated over frequency in the 
range between 10 Hz < / < 14 Hz). To facilitate compari- 
son with theory, the simulations are performed with periodic 
boundary conditions in the horizontal plane, otherwise the 
parameters are the same as in Fig. [3l 



Brillouin zone, where the intersection occurs (note also 
images of the hybrid mode outside the first zone). The 
theoretical curves calculated from Eqs ([I]), ([2]), (|A1|) . and 
()A2p for the same set of parameters demonstrate excel- 
lent agreement with simulations. 

The reason why the fluctuation intensity is anoma- 
lously high at the position of the hybrid mode is the 
heating induced by the wake-mediated coupling. Fig- 
ure [3] represents the marginally unstable regime, when 
Im cjhyb(— Im l^hyb " ^^) ^ 0- ^his regime (also ob- 
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served in experiment—) the effective growth rate of the 
instabihty is low, which ahows us to obtain the fluctu- 
ation spectra before the crystal is eventually destroyed. 
It is worth noting, however, that the hybrid mode would 
usually appear as a "hot spot" even if the instability is 
suppressed by friction and Im cjhyb < 0. This is because 
the heating due to mode coupling is present also in stable 
regime. It becomes negligible only when Im Qhyh ^ \^-> 
and then the fluctuation intensity coincides with the in- 
tensity of DL modes outside the hybrid zone (where it is 
determined by the background neutrals or temperature 
of thermostat in simulations). 

At the initial stage of the instability (and, of course, 
when the instability is suppressed) the effect of nonlin- 
earity on the DL modes is not significant. In this case 
the contour plot of the growth rate Im co'hyb(k) reason- 
ably reproduces the distribution of fluctuation intensity 
in the vicinity of the hybrid mode. This is illustrated in 
Fig. m where we plotted Im cjhyb(k) predicted by theory 
and the (frequency-integrated) fluctuation spectra from 
simulations. At a later stage, however, a nonlinear cou- 
pling between different modes becomes essential. This 
results in a variety of new phenomena, such as the en- 
ergy cascades from the "hot spots" in the k space and 
the generation of secondary harmonics at the double fre- 
quency of the hybrid mode (which are clearly seen in 
Fig. 3 of Ref.^^). Nevertheless, since the principal aim 
of this paper is to focus on distinct features of the linear 
mode coupling^ we leave aside the discussion of numerous 
nonlinear coupling effects. 

For the plots shown in Fig. |4]we used the same set of 
parameters as for Fig. [3l but for simulations we assumed 
periodic boundary conditions. The latter not only al- 
lowed us to directly compare the theory and simulations, 
but also to probe the effect of weak deviation from pe- 
riodicity caused by the horizontal confinement. We see 
that the only noticeable difference between Figs [3b and 
IDd is the intensity of the hybrid mode. This reflects the 
fact that the growth rate of the instability critically de- 
pends on the interparticle distance A (resulting in expo- 
nentially growing difference in the heating at the initial 
stage). Therefore in the case of parabolic confinement 
(where A in the center is somewhat smaller than the av- 
erage value) the intensity of the hybrid mode is higher 
than in the homogeneous case shown in Figs IDd. Oth- 
erwise, the positions and the overall shape of the "hot 
spots" in these two plots practically coincide, suggest- 
ing that weak inhomogeneity which is always present in 
experiments due to the horizontal confinement does not 
really affect the periodic structure of the Brillouin zones. 

Finally, to illustrate the difference between "shallow" 
and "deep" intersections of the out-of-plane and longitu- 
dinal in-plane modes, in Fig. [5] we demonstrate the fluc- 
tuation spectra obtained for the eigenfrequency of ver- 
tical confinement /^fonf which is slightly lower than that 
in Fig. [3l This naturally shifts the out-of-plane mode 
down and thus makes the resulting hybrid mode much 
more pronounced: Figure [5] shows that for deep intersec- 
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FIG. 5. Illustration of "deep" mode intersection (simula- 
tions). The eigenfrequency of the vertical confinement is 
fconi ~ 13.12 Hz (spectra are integrated over frequency in 
the range between 9 Hz < / < 13 Hz), the damping rate is 
ly = 2.19 s~^. Otherwise the parameters are the same as in 
Fig. El 

tions the separate "islands" of the hybrid mode merge 
into a joint broad "ring" (within the first Brillouin zone) 
and the "hot" region occupies substantially larger part 
of the k-plane. The magnitude of Im l^hyb for deep in- 
tersections is usually much larger than for the shallow 
intersections and therefore the threshold conditions to 
trigger the mode-coupling instability in the former case 
are substantially relaxed (see discussion in Sec. IIIIPp . 
Note that in this case, in order to keep the crystal in the 
marginally unstable regime (which allows us to recover 
the fluctuation spectra) we increased the damping rate 
toiy = 2.19 s-^ 



B. Angular dependence 

The dependence of DL modes on the direction of k 
is evident from the discussion above. For the hybridiza- 
tion and the onset of the mode-coupling instability this 
dependence becomes crucial and therefore it deserves a 
separate discussion. 

Figures [3] and m show that the first crossing of the out- 
of-plane and longitudinal in-plane modes always occurs 
dit 30°, at the border of the first Brillouin zone. The 
area of the hybrid mode in the k plane rapidly grows 
as the intersection gets deeper, and eventually the hy- 
brid mode exists for any 6>, as in Fig. [5l The cross- 
ing depth is controlled by changing the particle density 
or /and the eigenfrequency of the vertical confinement 
(see Sec Hiroi) . 
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Therefore, when the DL fluctuation spectra are mea- 
sured in experiments, under conditions close to the onset 
of hybridization (and, if not suppressed - of the mode- 
couphng instabihty), special attention should be paid to 
the orientation of k with respect to the crystalline lattice. 
Under these circumstances the direction should be chosen 
according to the purpose of the experiment: If detection 
of the mode coupling is the goal, then 9 = 30° is the best 
choice; if unperturbed modes are to be measured, then 
^ = 0° should be used. The experimental technique must 
therefore provide angular resolution, which is easily done 
for the in-plane modes, but is not always possible for the 
out-of-plane mode (since it is difficult to achieve angu- 
lar positioning of the plasma crystal in an experiment). 
For instance, if the out-of-plane mode is measured with 
a flxed side-view camera^^ then the angle is set during 
the experiment and cannot be adjusted at the analysis 
stage, but if it is measured with a top- view camera (from 
the variation of individual particle intensity)^ then 6- 
resolved wave spectra can be obtained. 



C. Mixed polarization 




FIG. 6. Mixed polarization near the hybrid mode (simula- 
tions). Shown are the fluctuation spectra with the transverse 
out-of-plane (a) and longitudinal (b) polarization, for k at 
= 30° and parameters of Fig. (3] Arrows indicate traces of 
mixed polarization. Far from the hybrid mode ("hot spots") 
the modes become purely out-of-plane (a) and in-plane (b). 

As we pointed out in Sec. |TT1 the DL wave modes 
are independent only when the wake-mediated coupling 
can be neglected: Far from the intersection or if the 
modes do not cross at all, the out-of-plane mode is 
purely transverse, one of the in-plane modes is longi- 
tudinal and the other is transverse. However, near the 
hybrid mode, where the coupling between the flrst two 
modes is strong, the theory predicts the emergence of 
mixed polarization^^: As the modes approach the junc- 



tion point, their polarizations become oblique (in the 
plane formed by k and the vertical axis), but the eigen- 
vectors remain mutually orthogonal. At the point where 
they merge the polarization exhibits a discontinuity - 
the eigenvectors of the two modes form a single (hybrid) 
eigenvector with the elliptic polarization. 

Traces of the mixed polarization for the out-of-plane 
mode have been observed in our recent experiments^ 
(supplemental Fig. S2 therein) as well as in the exper- 
iment by Liu et alS^ (shown in their Fig. 4b and er- 
roneously interpreted as a "new mode"). In Fig. [6] we 
illustrate the mixed polarization observed in simulations 
(performed for the conditions of our experimenti^i) . The 
magnitude of this effect (oc < 1) is quite small for typi- 
cal experimental conditions, so that its detection requires 
very careful analysis and high signal-to- noise ratio. 



D. Thresholds 

The mode-coupling theory predicts two distinct thresh- 
olds: (i) Confinement threshold^ which determines when 
the out-of-plane and longitudinal in-plane modes inter- 
sect; it depends on the combination of the vertical con- 
finement frequency and the lattice constant (particle 
number density). This threshold identifies the sufficient 
condition for the hybrid mode to form, but only the nec- 
essary condition for the instability to set in. The suffi- 
cient condition for the instability is determined by (ii) 
damping threshold - this identifies the critical value of 
the damping rate (controlled by pressure) below which 
the instability is not suppressed. Let us discuss thresh- 
olds (i) and (ii) separately (we recall that all frequencies 
are normalized by I^dl)- 



1. Confinement threshold 

The first intersection of the out-of-plane and longitu- 
dinal in-plane modes takes place at the border of the 
Brillouin zone^, for k at ^ = 30°. Then, using the ex- 
pressions for r^v(k) and Vt\^\\ (k) we readily obtain the con- 
dition for the hybrid mode to form (confinement thresh- 
old). 



k=ke 



(8) 



where |kcr| = ^7 



.(cr) 



IS 



Generally, the critical confinement frequency ^^^Q^f 
a rather comphcated function of tz^ 6, and q; it is given in 
Appendix [HI Eq. (jPip . However, Eq. (|Aip suggests that 
to the accuracy 0(6) the r.h.s. of Eq. (j8]) is independent 
of S and is proportional to — q. Thus, in this limit 

(cr) 



of "small" 5 the combination ^conWl ~ Q 



^^conf 1^ 

a function of k only (i.e., of the particle density) which 
is given by Eq. ()D2p . Figure [T^i shows this universal 
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FIG. 7. Universal mode-coupling thresholds (limit of small 
6). (a) Critical value of the vertical confinement frequency, 
^ionf (confinement threshold, upper curves), which identifies 
the formation of the hybrid mode, shown versus the screen- 
ing parameter k = A/ A. Also plotted is the minimum value 
^i^nT'* (below which only multilayer crystals can exist, lower 
curves), (b) Critical value of the damping rate, z^*^^^-* (damp- 
ing threshold), which determines the sufficient condition for 
the instability to set in, shown versus the confinement control 
parameter, 11 = (^conf)^ ~ ^conf? for k = 0.5, 1, 1.5, and 2. 
The upper bound for the critical damping rate, 

^(cr,max)^ is 

reached when Qconf = ^i^nT'* (marked by bullets). The inset 
in (b) shows ^y(^^'™^^) for different the mode-coupling in- 
stability is always suppressed above the dotted line. In both 
panels, solid lines are exact 2D results derived in this paper, 
dashed lines represent results of a ID NN theory^^. Tilde 
denotes universal normalization (see text for details). 
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FIG. 8. Role of "finite" S. Shown is the relative deviation 
of ^^^onf (solid line) and ^p^nT'* (dashed line) from the values 
calculated in the limit of small S (see Fig. [7^). The results 
are illustrated ioi k = 1 and three different values of q. 



dependence in the useful range 0.5 < < 2 (upper solid 
line). 

One should always keep in mind that l^conf cannot be 
arbitrarily small^ii^: The out-of-plane mode ^^v(^) = 
^conf ~ decreases monotonically with |k| and attains 
a minimum at the border of the first Brillouin zone. If 
< this formally implies instability. Physically, this is 
because the vertical confinement in this case becomes too 
weak to keep the particles in the monolayer and the bi- 
furcation to a multilayer crystal occurs (i.e., a monolayer 
crystal is no longer a ground state). Thus, irrespective 
of the value of damping rate, a monolayer plasma crys- 
tal can only exist if the vertical confinement frequency 
exceeds the minimum value defined by^ 

This value is given by Eqs (|D3|) and (|D4|) for arbitrary 
and small values of (5, respectively. In the latter case we 

obtain a universal dependence ~ Q = ^il^nf^ 

on which is also plotted in Fig. [7^ (lower solid line). 

The effect of "finite" S [i.e., beyond the accuracy 0{S)] 
is illustrated in Fig. [8] for = 1 (results are practically 
independent of tz in the range considered here). We plot 

the relative deviations of ^^^onf ^conf ^ from the values 
given by Eqs (|D2p and (|D4p , respectively. The deviations 
increase as oc for ^ < 0.3 and then rapidly saturate at 
the level (1 -g)-^/^ - 1. [The latter is because both 
and ^coi!f depend on the wake parameters for 
S ^ 1^ whereas at <C 1 they scale as oc ^/l — q.] Note 
that the curves have almost identical shape at different q^ 
which indicates that the shown correction approximately 
scales as oc (1 — ^1 — q). Given typical experimental 
errors (see Table |I] and Fig. [JO]), we conclude that the 
corrections shown in Fig. [8] can generally be neglected. 
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i.e., the accuracy of universal dependencies obtained in 
the hmit of smah S is sufficient for practical purposes. In 
some cases, however, the effect of "finite" 6 can never- 
theless play an important role (see Sec. IIVB]) . 

Finally, let us compare the "exact" 2D results (in the 
limit of small S) with the results of approximate ID 
NN mode-coupling theory summarized in Sec. Ill CI For 
the ID case, the critical and minimum confinement fre- 
quencies are derived from Eq. ([5j) using the conditions 
^h(^) = = respectively. For the 

comparison with the 2D results, we multiply the ID re- 
sults by constant "form factors" reflecting mismatch of 
ID and 2D environments, which yields 



conf 

(min) 
conf 



2F 



2F^ 



^T^y(^ 

(cr) 



3/^ 



3/^ ^)e 



(10) 



^)e- 



The form factors F^"^^^^ = 1.25 and f'^^^^^ = 1.49 are 



determined from a fit to the 2D results, Eqs (|D2[) and 
(|D4p . The dependencies given by Eq. (p!Q|) are shown in 
Fig. [TK by dashed lines. 

Thus, we see that the ID NN theory provides excellent 
qualitative description of the confinement thresholds - 
even for small a^; by multiplying these results with con- 
stant factors about unity we also get very good quanti- 
tative agreement with 2D results. 



2. Damping threshold 



As pointed out in Sec. IIIB| for = 30° (as well as 
for 0°) the transverse in-plane mode becomes exactly de- 
coupled. In this case the out-of-plane and longitudinal 
in-plane modes are governed by the following equation 
[similar to that for the ID NN model, Eq. ^]\ 

In the vicinity of the intersection (l^o^ko) the coupled 
modes are described by Eq. (|4]) with e :^ (cry/2r^o)^- 
The group velocities Uhy and Uv are collinear with k 
due to the azimuthal symmetry at 6> = 30°, so that the 
hybrid mode is given by 

l]hyb(k) = 1^0 + ^(Uhll + Uv) • (k - ko) 

1 



±-zy/a2(ko)/l^g - |Uh|| - Uv|2|k - koP. 

Thus, the width of the hybrid mode is |k — ko| = 
cr^y (ko )/ (1^0 1 Uh II — Uv|)~"^, and the maximum of 
Im l]hyb(k) is attained at the crossing point ko. From 
this we immediately obtain the sufficient condition for 
the instability (damping threshold), 

(11) 

where VLq and ko are determined from the crossing con- 
dition. 



^<z/(^^) =a^(ko)/l^o, 



^0 = vohT^ 



k=ko 



2 

conf 



2a^ 



(12) 



In order to relate Eq. ([TT]) to the necessary instability 
condition [confinement threshold, Eq. (j8j)] we introduce 
the confinement control parameter^ 



n = (^^i:„\)^ 



(cr) x2 



which varies in the range from zero to (^conf) 
(^lorif ^)^- is obtained from Eq. ([Hj), 



n 



(ah + 2av + /3)|p 



(13) 



where the r.h.s. is calculated as the difference at kcr and 
ko. Thus, Eqs (pTj) and (p!3|) determine the dependence of 
the critical damping rate i^^^^^ on the control parameter 
II via the parametric dependence on |ko|. For arbitrary 
(5, this dependence is derived in Appendix [D) Eqs (|D5[) 
and jnH). 

Similar to the confinement threshold discussed in the 
previous subsection, for the damping threshold one can 
also obtain a universal dependence in the limit of small 5. 
It is calculated in Eqs (|D7|) and (|D8[) . where the r.h.s. are 
functions of n only (along with the parameter |ko|) and 
the dependence on the wake parameters is explicitly given 
by the l.h.s. By introducing the universal normalization 
for the critical damping rate and the confinement control 
parameter. 



qS 



_i^(cr) 



n = 



n 



1 



(14) 



we conclude that z>^*^^^ is a universal function of k, and IT. 
Figure [Jb shows z>^^^^ versus IT for four characteristic val- 
ues of hi (solid lines). For "shallow" intersections, when 
^conf is slightly below the critical value ^conf (^) 

ical damping rate increases as z>^^^^ oc a/iI, which can be 
readily derived from Eqs (|D7|) and (|D8|) . Thus, for given 



^(cr) 



and n (with ^conf ^ ^ ^conf < ^conf ) monolayer plasma 
crystals are unstable if < {qS/y/1 — q)j)^^'^^ (/^, 11) and 
stable otherwise. 

The critical damping rate is always bound from above, 
< z>(cr,max) a givcu A^, the upper bound 

~(cr,max)j^^j COrrCSpOUds tO l^conf = ^1™^^) ^^^^ 

marked by bullets in Fig. [Tb. These points shown 
in the inset as a function of n identify the absolutely 
stable region: The mode-coupling instability in mono- 
layer crystals is always suppressed when v exceeds 

(gVvT^)^^^''''"'''''^^)- 

According to the ID NN theory [see Appendix iB) Eq. 

(|B1[) ] the critical damping rate is 



j>(cr) 



3/^ - 



^2 + 2/^ + 2 



i.e., it also scales as a square root of the control parame- 
ter. The upper bound for i^^^^^(/^) is obtained by substi- 
tuting Eq. (fTQ|) for the maximum of 11, which yields 



k=ko 



~(cr,max) _ 2^^^^'"^^^ V(^"^ + ^A^-^ + e"'^, (15) 
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TABLE I. Parameters of 2D crystals for different experiments (marked I to XII). The mode crossing is indicated as observed 
in that study. 
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The value is too high for the mode crossing (see Fig. [TQ|) . amended in Erratum^. 



* For these data, error bars are available: ±15 % for k and A, ±10 % for |Q|, ±8 % for f^ll^, and ±4 % for A. 



where the corresponding form factor is i^(cr,max) _ 28. 
Equation ([IS]), which is plotted in the inset of Fig. [Jt) by 
the dotted hne, also provides excellent agreement with 
the 2D model. Hence, this is a convenient formula to 
identify conditions of absolute stability for a 2D crystal. 



IV. COMPARISON WITH EXPERIMENTS 

Wake-induced mode coupling has been observed in sev- 
eral experiment^— However, the analysis of exper- 
imental results can be fastidious and needs to be done 
carefully. The fingerprints of the mode-coupling must 
be identified properly in order to avoid misinterpreta- 
tions of the results. Therefore, in this Section we use 
the experiment of Ref. as a characteristic example to 
demonstrate the principal features of the wake-induced 
mode coupling. Then, we give an overview of literature 
currently available on experiments with 2D crystals (sta- 
ble and unstable)^ —i^i^, the summary of experimen- 
tal parameters is given in Table HI Finally, we discuss the 
stability of crystals in terms of thresholds introduced in 
the previous Section. 

One of the most clear experimental evidence of mode- 
coupling in a 2D plasma crystal was reported in Ref^^ 
(experiment H in Table H]). A monolayer of microparticles 
was formed by levitating them in a plasma sheath above 
the lower rf electrode in a capacitively coupled discharge 
in a (modified) GEC chamber. The microspheres were 
made of melamine-formaldehyde and had a diameter of 
9.19 ± 0.09 fim. The low argon pressure of 0.76 Pa en- 
sured that the particle motion was only slightly damped. 
A high quality of the monolayer (with no detected par- 
ticles above or below) was verified by using a side-view 
camera. Details of the experimental procedure are given 



in Refs^^. 

The experimentally measured fluctuation spectra of 
two in-plane DL modes are shown in Fig.[9l Since this ex- 
periment was performed at conditions close to the onset 
of instability, all essential fingerprints of mode coupling 
are clearly visible: 

• The mode coupling critically depends on the direc- 
tion of the wave vector k. The comparison of wave 
spectra calculated for l9 = 0° (Fig. ^) and = 30° 
(Fig. [SJd) reveals that the onset of mode coupling 
first occurs at 6> = 30° - this results in the promi- 
nent "hot spots" in Fig. [SJd. Colored red, they are 
almost two orders of magnitude more intense than 
the "normal" branches of the fluctuation spectra. 
The hot spot that lies within the first Brillouin zone 
(left one) represents the hybrid mode, the other one 
is its image in the higher Brillouin zone. These ob- 
servations are in agreement with theory and com- 
plimentary simulations (Fig. [3]), which predicts that 
the mode coupling first takes place at 6> = 30° for 
"shallow" mode intersection. 

• The less intense spot seen in Fig. [9^ at 6> = 0° 
and kA 6.5 is outside the first Brillouin zone (in 
Fig. [Hh the location is along /c^^-axis, at the right 
edge). As explained in Sec. IHDi this spot should 
be mapped into the first Brillouin zone, where it 
represents the gap between two hot spots of the 
hybrid mode (in Fig. [Sh, along the border of the 
Brillouin zone at :^ 30°), and therefore its inten- 
sity is almost one order of magnitude smaller than 
the intensity of the hot spots. 

• Traces of mixed polarization (out-of-plane mode) 
are clearly seen next to the hybrid mode in Figs. [9^ 
and[9]3. This gives a direct experimental proof that 
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FIG. 9. In-plane fluctuation spectra for "shallow" mode intersection (experiment IJ^ in Table |T}. Shown are the DL modes 
(positive branches) for k at (a) ^ = 0° and (b) = 30° , and (c) spectrum in the k-plane integrated over frequency (in the range 
between 10.5 Hz < f < 14.5 Hz). The dashed lines show the border of the first Brillouin zones, solid lines in (a) and (b) are 
theoretical curves for all tree principal DL modes (assuming q — 0.3 and S = 0.33). Fluctuation spectra in (a) and (b) reveal 
clear traced of the out-of-plane mode (mixed polarization). In order to reduce noise in (c), the spectrum is an average of six 
spectra corresponding to equivalent directions in the crystal (using the invariance of hexagonal lattice to 60° rotations). 



the out-of-plane mode becomes oblique close to the 
intersection point. A discrepancy of 5 % between 
the experimentally measured in-plane projection of 
the out-of-plane mode and its theoretical dispersion 
relation can be explained by "finite-^" correction to 
the theory shown in Fig. [8] (which is about 4 % for 
the chosen values of wake parameters). 

• Finally, all "hot spots" predicted by the theory 
and observed in simulations can be simultaneously 
seen in the experimental wave spectrum in k-plane, 
which is shown in Fig. [St. 

Next, we analyze the stability of 2D plasma crystals 
in terms of the confinement threshold (which determines 
when the hybrid mode emerges). Available literature 
data on stable and unstable crystals is summarized in Ta- 
ble m and the measured vertical confinement frequency 

^'^/ionf compared to the theoretical threshold ^^^onf 
Fig. [TOl Note that for the comparison we change the 
normalization of 27r/^^^^£: Instead of using ftj^L - a nat- 
ural frequency scale in theory, now we employ the 2D 



analogue of dust plasma frequency. 

This normalization allows us to minimize the experi- 
mental errors, since the interparticle distance is usually 
measured with much higher accuracy than the screening 
length (see footnote to Table H]). 

Figure [To] shows that all experiments listed in Table |I] as 
having no mode crossing are indeed represented by points 
well above the confinement threshold line. To the con- 
trary, the points for experiments marked by mode cross- 
ing are located very close to the threshold line or below it. 
The only exception is the experiment VIII^^ by Liu et al: 
The wake-induced mode coupling was clearly seen in the 
fluctuation spectra measured in^^, but the reported ex- 
perimental parameters (yielding '^^f^J^f/^o — 6, marked 
by green up triangle in Fig. fTO]) certainly do not allow the 
mode crossing, which indicates large measurement errors 
(amended in the subsequent Erratum^^). This example 
demonstrates that, by comparing experimental param- 
eters to the theoretical confinement threshold near the 



12 





I'll 


1 

• ▲ 




10 






No mode 


8 






crossing 


o 

S 


1 

T 


• 




^ o 

M— 

4 


• 

' ~— — ^ ^ 










Hybrid 










mode 


2 


Multilayer 

1 1 1 1 


1 





0.5 1 1.5 2 



K 

FIG. 10. Comparison of the vertical confinement frequency 
measured in experiments (Table |l|) with the theoretical con- 
finement threshold (Fig. [7^): Blue diamonds are for experi- 
ments I and II, red circles for III- VI, green up triangles for 
VII and VIII, open green up triangle for IX, blue square for 
X, and purple down triangles for XI and XII. The vertical 
resonance frequency is normalized to ujq — \J Q'^ /m/S.^ . For 
experiments I to VI, error bars are available: ±15% for ^t, 
±8% for f'^ll^, and ±12% for cjq. 

point where the hybrid mode (hot spot) emerges in the 
fluctuation spectra one can perform very sensitive con- 
sistency test for the deduced parameters. 



A. Effective charge and wake parameters 

At this point, let us briefly discuss the relation between 
the actual charge carried by particles and the effective 
charge which determines the dispersion properties of DL 
waves in 2D plasma crystals. 

The common method to determine Q and A in 2D crys- 
tals is to fit the measured fluctuation spectra for the in- 
plane modes by theoretical curves. The method is based 
on the assumption that the interparticle interaction po- 
tential can be approximated by the Yukawa form3 and 
provides rather high accuracy. However, the dispersion 
relations are determined by the combined effect of direct 
interparticle interactions and the wake-mediated inter- 
actions. While the former contribution scales as oc Q^, 
the latter depends on a particular form of the wake. For 
the model of a point-like wake charge used in this pa- 
per the wake- mediated interactions scale as ex gQ, which 
modifies the dispersion elements of the dynamical matrix 
(|A1[) . In the limit of small 5 (which turns out to be a re- 
markably accurate approximation, see Fig. [8]) this yields 
the common renormalizing factor (1 — q) in Eq. ()A1[) . 
and therefore the dispersion relations are determined by 
the effective charge Qeff = — Q Q- 



This fact should always be kept in mind when fit- 
ting the fluctuation spectra by theoretical dispersion re- 
lations: In theory, frequencies are normalized by I^dl oc 
\Q\ and hence the explicit dependence on the factor 
{1 — q) naturally disappears. This allows us, e.g., to 
directly compare the experimentally measured values of 
^'^/conf/^o with the theoretical confinement threshold 
without knowing the wake parameters (see Fig. [TQ|). 

Of course, the comparison of theory and experiments 
in terms of the damping threshold requires knowledge 
of q and 6, since the growth rate of the hybrid mode 
is directly proportional to the effective "dipole moment" 
of wake qS. As was shown in Sec. IIII D 2[ in order to 
compare the measured value of u/ujo (last column of Ta- 
ble H]) with the theoretical damping threshold plotted in 
Fig. [TJd, we need to known the factor qS/{l — q) [see Eq. 
(dl]); additional factor — q in the denominator here 
is because ujq in experiments is determined by Qeff]- In 
experiments, however, the wake parameters are generally 
unknown, so that a reasonable assumption should always 
be made in order to perform the analysis in terms of the 
damping threshold. 

B. Effect of "finite" ~5 

As we mentioned in Sec. IIIID 2\ under certain con- 
ditions the "finite-(5" correction shown in Fig. [8] can be- 
come important. In particular, this happens in "marginal 
cases" , when the normalized value of the confinement fre- 
quency is very close to the corresponding threshold line 
(viz., when the distance to the line is comparable or less 
than the correction). For instance, such situation might 
occur for experiments II, VI, IX, and X shown in Fig. [TOl 
However, typical experimental errors are relatively large, 
so that accounting for the corrections is hardly required 
in these cases. 

As regards the comparison of theory with numerical 
simulations, where all input parameters are known pre- 
cisely, the "finite-(5" corrections can become crucial. To 
illustrate this, let us consider simulations of the shallow 
mode intersection shown in Fig. [3l The (normalized) 
confinement frequency in this case is ~ 1 % larger than 
^conf according to the universal criterion for the 
confinement threshold the system is still marginally sta- 
ble. However, by taking into account the correction to 
the confinement threshold from Fig. [8] (which is about 
±4 % in this case) we conclude that the system is in fact 
marginally unstable - as observed in the simulations. In 
contrast, for the case of the deep mode intersection shown 
in Fig.[5](where the confinement frequency is just 6.5 % 
smaller) the correction plays practically no role. 

V. CONCLUSIONS AND OUTLOOK 

The theory of mode-coupling instabilit y^^ provides 
detailed picture of a plasma- specific melting scenario op- 
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erating in 2D plasma crystals. The melting associated 
with the wake- mediated coupling between the longitu- 
dinal in-plane and out-of-plane modes can only be trig- 
gered if (i) the modes intersect and (ii) the neutral gas 
damping is sufficiently low. Even if the instability is sup- 
pressed by the damping and the melting does not occur, 
the coupling always results in the formation of a hybrid 
mode which is revealed in fluctuation spectra by anoma- 
lous "hot spots" emerging at distinct positions in the 
k-plane. In the vicinity of the hybrid mode, one can ob- 
serve traces of mixed polarization for the two intersecting 
wave branches. 

In this paper we showed that all these features can be 
considered as distinct fingerprints to identify the onset 
of the wake-induced mode coupling. This, in turn, al- 
lows us to determine certain combinations of the crystal 
parameters with a fairly good accuracy. For instance, 
the theory predicts a well-defined confinement threshold 
(for the ratio of the vertical confinement frequency to the 
dust-lattice frequency) at which the "hot spots" emerge. 
It is noteworthy that the effect of wakes on the disper- 
sion relations is automatically taken into account via the 
charge renormalization, which yields the effective charge 
- the value which is deduced from the analysis of exper- 
imental fluctuation spectra. 

On the other hand, the damping threshold as well as 
the evolution of the mode-coupling instability (if not sup- 
pressed) are determined by parameters of the wake - 
which are still poorly known for experiments. By identi- 
fying the onset of the instability and comparing it with 
the damping threshold derived above, one could obtain 
effective "dipole moment" of the wake. Furthermore, one 
could measure the energy growth rate at the initial in- 
stability stage and compare it with the theoretical pre- 
diction. Systematic studies in this direction would be 
highly desirable. However, we note that such measure- 
ments would only give us some mean characteristics of 
the wake, e.g., the first ("dipole") moment in the multi- 
pole expansion. 

We believe that future research into the stability of 
2D plasma crystals should be focused, in particular, on 
the implementation of self-consistent wake mode b^^i^^i^^ 
for the calculation of the mode coupling (which should 
also account for the spatial variation of the interaction 
parameters'*^) . Such studies will allow us to gain 
deeper insight into the plasma-specific mechanisms of 
melting operating in 2D systems and provide more reli- 
able basis for future research into generic melting mech- 
anisms. 



following auxiliary functions: 
^(x) = {x~^ + x~ 



+ x-3)e-^. 



S(x) {x-^ + 3x-2 + 3x-^)e- 
A(x) = i[S(x)-^(x)]. 



Then the dispersion elements of Eq. ([T]) are given by the 
following sums over integer m and n with excluded (0, 0): 

m,n ^ '- ^ 

X sin^ • s, 

m,n ^ '- ^ 

X sin^ • s, 

m,n 

X {Sy - sl)s~'^ sin^ ^k • s. 



sin^ ^k • s. 



(Al) 

Here, the vector s with the components Sx = hm-\-n and 



Sy = and the absolute value s 



\Jvn? + mn + v? 
characterizes the relative positions of all neighbors in the 

hexagonal lattice; we also introduced = Vs"^ 
The coupling elements are given by 

o-x,y = ^{i^ss)sx,ySj'^ sink • s. (A2) 



In the NN approximation (formally applicable for ^ 1), 
the modes ^^v(^) ^h||(^) (whose coupling causes the 
hybridization) obtained from Eqs ()Aip and ()A2p coincide 
with Eqs. (j5j)-([7]) for a ID string. In the limits q ^ or 
^ ^ oo, Eqs ()Aip and ()A2p reduce to conventional (wake- 
free) results for DL modes in 2D crystals. For S ^ and 
finite q one also obtains the wake- free results, but with 
the particle charge renormalized by the factor ^/l — q. 



Appendix B: Hybrid mode for ID model 

The intersection point between the horizontal and ver- 
tical modes, (l^c^o), is readily derived from Eq. (j5j): 



Appendix A: Elements of dynamical matrix D [Eq. ([T]l] 

In Ref^^ we calculated the elements of D with the accu- 
racy 0{5). Although usually 5 is small indeed, such a lin- 
ear expansion might not be always sufficient for some ex- 
perimental conditions (see Sec. IIVB]) . Below we present 
the results for arbitrary 5. For brevity, we introduce the 



I tl^ -\- 2k -\- 2 .1 ^conf 

^^0 = \ 9 , o — ^^conf, sm -ko 



k'^ -\- 3k -\- 3k 



(cr) ' 
conf 



where ^^onf given by the first Eq. (fTQ|) with i^^^onf ~ 
1. The maximum growth rate of the hybrid mode is 
obtained from Eq. (j6]), 

a{ko) 



Im r^hyb(^o) 



21^0 



(Bl) 
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where a{k) is given by Eq. ([7]). 
damping rate z/^'^^^ = <j(/co)/^o- 



This yields the critical is 



■ exp 



Appendix C: Molecular dynamics simulations 



The charged particles were confined in a vertical 
parabolic well with an eigenfrequency which was varied 
in the range of 15 Hz < f^^l^f < 25 Hz. In the horizontal 
plane, two options were investigated: parabolic confine- 
ment with an eigenfrequency f^^^^ < 0.5 Hz and periodic 



boundary conditions (with f^^^^ 
motion for each particle is 



0). The equation of 



where Vij is the distance from the particle i to the particle 
j, T^.. is the distance to the wake of the particle j, Q is 
the particle charge, and q is the point-like charge of the 
wake placed at a distance S below each particle. The 
Langevin force is defined as 

{U{t)) = 0, (L,(t + r)L,(t)) = 2iymTS^jS{r), 

where T is the temperature of the thermostat (here, S is 
the Kronecker delta or the Dirac delta function). 

In order to simulate systems of a size comparable to ex- 
perimental (with TV > 10^ particles) without any approx- 
imation (i.e., without cutoff radius in the interactions), 
the code was written in CUDA C and run on NVIDIA 
Tesla C1060 GPU or C2050 GPU computing cards. The 
equations of motion were solved using Beaman algorithm 
with predictor-corrector. 



Appendix D: Calculation of thresholds for 2D model 



where m is the particle mass, u is the damping rate, 
the position of the ith particle, and C(r^) represents the 
confinement force. The force of interparticle interaction 

I 



To derive the confinement threshold, we substitute the 
dispersion elements from Eq. (|Aip in Eq. (|8|) which 
yields 



(cr) 
conf 



I'^Vm |(1 +P) [S(/^5) -qE{K.ss)] ^qPs^'^ {3^p)E{K.ss)y 



(Dl) 



r 



Here the fact that |kcr|5y = uitt is taken into account, 
so that it is convenient to introduce the parameter r]rn = 
^ [1 — (—1)^] (which is equal to 1 or for odd or even 
m, respectively). Furthermore, for brevity we use p = 



For the accuracy 0((5), we can neglect in Eq. (jPip the 
second term in curly braces and set ss = which results 

I 



in the following expression: 



(cr) 
conf 



,??m(l+p)5(Ks), 



(D2) 



where the r.h.s. is a function of only. Thus, Eq. (|D2| 

determine the universal dependence of ^^^onf ^• 

The minimum value of the confinement frequency de- 
fined in Eq. (j9j) is 



(D3) 



for the accuracy 0{S) it is reduced to 



(min) 



7r^ 



(D4) 



For the damping threshold, we substitute the coupling 
element [Eq. (|A2|) ] in Eq. (pTj) and the dispersion ele- 
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ments in Eqs (p!2]) and ([13]) . After some algebra we get 



= ^ , (D5) 

5] ] *(ks) - + P [S(ks) - q^{KSs)] + qd^Sg ^ {I + p) ^{kss) \ sin^ ^\^o\Sy 



n = ^ (^m - sin^ WM^y) {(1 +P) [^{^s) - q^{^ss)] + qS^sf (3 + p) ^(a^s^)} . (D6) 



I 

For the accuracy 0((5), we neglect the last terms in curly 
braces and set S5 = s. Then Eqs ()D5p and ()D6p are 
reduced to 



VI - (cr) 



^ s '^Sy'E{K.s)sm\]<iQ\sy 

^(crj ^ =^ (D7) 



^ (r?™ - sin^ \\ko\Sy) [l+p] 5(ks), (D8) 



n 



where the r.h.s. depend only on n and |ko|. Thus, Eqs 
(|D7p and (|D8p determine the universal dependence of 
the critical damping rate z^^^^^ on the confinement control 
parameter 11 via the parametric dependence on |ko|. 
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